% By Martin M. Andreasen, April 22 2010
% This function simulates the model when solved up to third order. The purning scheme is used.
% These codes are faster for a third order approximation because we use the kronecker
% product to eleminate loops
function [Y_sim,yields_sim,ER_sim,X_sim] = sim3rdPrun(model,shocks,orderApp)

% Unfold model
gx  = model.gx;
%gxx = model.gxx;
%gxxx= model.gxxx;
gss = model.gss;
gssx= model.gssx;
GGxxtil  = model.GGxxtil;
GGxxxtil = model.GGxxxtil;

hx  = model.hx;
%hxx = model.hxx;
%hxxx= model.hxxx;
hss = model.hss;
hssx= model.hssx;
HHxxtil  = model.HHxxtil;
HHxxxtil = model.HHxxxtil;

eta = model.eta;
nx  = model.nx;
sig = 1;

numSim = size(shocks,2);    %Number of simulations

% Allocating memory
nMat       = length(model.by0);
Y_sim      = zeros(model.ny,numSim);
yields_sim = zeros(nMat,numSim);
ER_sim     = zeros(length(model.matConShortRate),numSim);
X_sim      = zeros(model.nx*orderApp,numSim);

% The simulation
xf    = zeros(nx,1);    %The first order terms
%xs    = zeros(nx,1);    %The second order terms
vec_xfxf  = (eye(model.nx^2)-kron(model.hx,model.hx))\reshape((model.eta*model.eta'),model.nx^2,1);
xs        = (eye(model.nx)-model.hx)\(model.HHxxtil*vec_xfxf+0.5*model.hss);
xrd   = zeros(nx,1);    %The third order terms

% Defining matrices
if orderApp == 1
    for t=1:numSim
        X_sim(1:nx,t)   = xf;
        Y_sim(:,t)      = gx*xf;
        yields_sim(:,t) = model.byx*xf;
        ER_sim(:,t)     = model.rx*xf;
        
        % The states next period
        xf_p  = hx*xf + eta*shocks(:,t);
        
        % Updating xf
        xf  = xf_p;
    end
elseif orderApp == 2
    for t=1:numSim
        X_sim(1:nx,t)      = xf;
        X_sim(nx+1:2*nx,t) = xs;
        AA                 = kron3(xf,xf);
        Y_sim(:,t)         = gx*(xf+xs) + GGxxtil*AA + 0.5*gss;
        yields_sim(:,t)    = model.byx*(xf+xs) + 0.5*model.byxx*AA + 0.5*model.byss;
        ER_sim(:,t)        = model.rx*(xf+xs)  + 0.5*model.rxx*AA  + 0.5*model.rss;
        
        % The states next period
        xf_p  = hx*xf + eta*shocks(:,t);
        xs_p  = hx*xs + HHxxtil*AA + 0.5*hss;
        
        % Updating xf, xs
        xf  = xf_p;
        xs  = xs_p;
    end
elseif orderApp == 3
    for t=1:numSim
        X_sim(1:nx,t)        = xf;
        X_sim(nx+1:2*nx,t)   = xs;
        X_sim(2*nx+1:3*nx,t) = xrd;
        AA              = kron3(xf,xf);
        AAA             = kron3(xf,AA);
        BB              = kron3(xf,xs);
        Y_sim(:,t)      = gx*(xf+xs+xrd) + GGxxtil*(AA+2*BB) + GGxxxtil*AAA + ...
            0.5*gss + 3/6*gssx*xf;
        yields_sim(:,t) = model.byx*(xf+xs+xrd) + 0.5*model.byxx*(AA+2*BB) + 1/6*model.byxxx*AAA + ...
            0.5*model.byss + 3/6*model.byssx*xf;
        ER_sim(:,t)    = model.rx*(xf+xs+xrd) + 0.5*model.rxx*(AA+2*BB)    + 1/6*model.rxxx*AAA  + ...
            0.5*model.rss  + 3/6*model.rssx*xf;
        
        % The states next period
        xf_p  = hx*xf + sig*eta*shocks(:,t);
        xs_p  = hx*xs + HHxxtil*AA + 0.5*hss;
        xrd_p = hx*xrd + 2*HHxxtil*BB + HHxxxtil*AAA + 3/6*hssx*xf;
        
        % Updating xf, xs, xrd
        xf  = xf_p;
        xs  = xs_p;
        xrd = xrd_p;
    end
end
end